##########################################################################################

library(data.table)
library(optparse)
library(Seurat)
library(dplyr)

##########################################################################################
option_list <- list(
    make_option(c("--rna_file"), type = "character"),
    make_option(c("--out_path"), type = "character") 
)

if(1!=1){

    ## 单细胞表达文件
    rna_file <- "~/20231121_singleMuti/results/qc_atac/testis_combined.annotationCellType.qc.Rdata"

    ## 输出
    out_path <- "~/20231121_singleMuti/results/qc_atac/"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

rna_file <- opt$rna_file
out_path <- opt$out_path

dir.create(out_path , recursive = T)

###########################################################################################
## 导入数据
a <- load(rna_file)

###########################################################################################
#########################################
## 构造表达矩阵，原始的，如果直接提取需要10h左右
## 下面步骤可以加速到半小时以内，将原始的data转化为数值型变量
exp_numeric <- apply(scrnat@assays$RNA@data , 1 , as.numeric)
rownames(exp_numeric) <- colnames(scrnat@assays$RNA@data)
exp_numeric <- t(exp_numeric)

sco_exp <- sapply(unique(scrnat$cell_type),function(x){
    print(x)
    sapply(unique(rownames(scrnat)),function(y){
        mean(as.numeric(as.vector(exp_numeric[y,which(scrnat$cell_type==x)])))
    })
})

sco_exp2 <- data.frame(sco_exp)
sco_exp2$gene <- rownames(sco_exp2)
sco_exp2 <- sco_exp2[,c( ncol(sco_exp2) , (1:ncol(sco_exp2)-1) )]

out_file <- paste0( out_path , "/GeneExpression.MeanByCellType.tsv" )
write.table( sco_exp2 , out_file , row.names = F , sep = "\t" )

#########################################
## 构造表达矩阵,填补后的,已经是数值型
sco_exp <- sapply(unique(scrnat$cell_type),function(x){
    print(x)
    sapply(unique(rownames(scrnat)),function(y){
        mean(as.numeric(as.vector(scrnat@assays$MAGIC_RNA@data[y,which(scrnat$cell_type==x)])))
    })
})

sco_exp2 <- data.frame(sco_exp)
sco_exp2$gene <- rownames(sco_exp2)
sco_exp2 <- sco_exp2[,c( ncol(sco_exp2) , (1:ncol(sco_exp2)-1) )]

out_file <- paste0( out_path , "/GeneExpression.MeanByCellType.magic.tsv" )
write.table( sco_exp2 , out_file , row.names = F , sep = "\t" )

#########################################
## 细胞分成三类
SSC_SPG <- c("SSC" , "Differenting&Differented SPG")
SPC <- c("Leptotene" , "Zygotene", "Patchytene",
    "Diplotene" , "Early stage of spermatids")
SPT <- c("Round&ElongateS.tids" , "Sperm")

scrnat$cell_type <- ifelse(scrnat$cell_type %in% SSC_SPG , "SSC_SPG" , scrnat$cell_type)
scrnat$cell_type <- ifelse(scrnat$cell_type %in% SPC , "SPC" , scrnat$cell_type)
scrnat$cell_type <- ifelse(scrnat$cell_type %in% SPT , "SPT" , scrnat$cell_type)

sco_exp <- sapply(unique(scrnat$cell_type),function(x){
    print(x)
    sapply(unique(rownames(scrnat)),function(y){
        mean(as.numeric(as.vector(exp_numeric[y,which(scrnat$cell_type==x)])))
    })
})

sco_exp2 <- data.frame(sco_exp)
sco_exp2$gene <- rownames(sco_exp2)
sco_exp2 <- sco_exp2[,c( ncol(sco_exp2) , (1:ncol(sco_exp2)-1) )]

out_file <- paste0( out_path , "/GeneExpression.MeanByCellType.Combine.tsv" )
write.table( sco_exp2 , out_file , row.names = F , sep = "\t" )
